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We present a dual geometrical worm algorithm for two-dimensional Ising models. The existence 
of such dual algorithms was first pointed out by Prokof 'ev and Svistunov^ . The algorithm is defined 
on the dual lattice and is formulated in terms of bond-variables and can therefore be generalized 
to other two-dimensional models that can be formulated in terms of bond- variables. We also dis- 
cuss two related algorithms formulated on the direct lattice, applicable in any dimension. These 
latter algorithms turn out to be less efficient but of considerable intrinsic interest. We show how 
such algorithms quite generally can be "directed" by minimizing the probability for the worms to 
erase themselves. Explicit proofs of detailed balance are given for all the algorithms. In terms of 
computational efficiency the dual geometrical worm algorithm is comparable to well known cluster 
algorithms such as the Swendsen-Wang and Wolff algorithms, however, it is quite different in struc- 
ture and allows for a very simple and efficient implementation. The dual algorithm also allows for 
a very elegant way of calculating the domain wall free energy. 

PACS numbers: 05.10.Ln, 05.20.-y, O2.70.Tt 



I. INTRODUCTION 



Over the recent decades many powerful Monte Carlo 
(MC) algorithms have been developed, greatly enhancing 
the scope and applicability of Monte Carlo techniques. 
In fact, it is quite likely that these algorithmic advances 
have, and will continue to have, a far greater impact on 
the predictive power of Monte Carlo simulations than 
advances in raw computational capacity, a point that is 
often overlooked. The continued development of such 
advanced algorithms is therefore very important. Here 
we shall mainly be concerned with MC algorithms suit- 
able for the study of lattice models described by classical 
statistical mechanics. Some of the most notable devel- 
opments in this field have been the development of clus- 
ter algorithms by Swendsen and WangS^ and by WolfF^. 
More recent developments include invaded cluster algo- 
rithmsi that self-adjust to the critical temperature, flat 
histogram methods^, focusing on the density of states 
and techniques performing Monte Carlo sampling of the 
high temperature series expansion of the partition func- 
tiorii using worm algorithms^. Two of us recently pro- 
posed a very efficient geometrical worm algorithm.^'* for 
the bosonic Hubbard model. In this algorithm variables 
are not updated at random but instead a "worm" is prop- 
agated through the lattice, at each step choosing a new 
site to visit with a probability proportional to the rela- 
tive (among the different sites considered) probability of 
changing the corresponding variable. The high efficiency 
of the algorithm stems from the fact that the worm is al- 
ways moved and it is only through the actual movement 
that the local environment - the local "geometry" - comes 
into play. The ideas underlying this algorithm are quite 



generally applicable and the present paper is concerned 
with their generalization to classical statistical mechanics 
models. 

As the canonical testing ground for new algorithms we 
consider the standard ferromagnetic Ising model in two 
dimensions defined by: 
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Here < i,j > denote the summation over nearest neigh- 
bor spins. It is well known that the critical temperature 
for this model in two dimensions is fc^Tc = 2J/log(l -I- 
V2) = J2. 26918 . . .. When investigating magnetic ma- 
terials modeled by classical statistical mechanics such as 
Eq. iQJ using Monte Carlo methods one has to take into 
account the effects of the non-zero auto-correlation time 
T (defined below) that is always present in Monte Carlo 
simulations. The auto-correlation time describes the cor- 
relation between observations of an observable O{to) and 
0(^0 + t Monte Carlo Sweeps (MCS) apart. The auto- 
correlation time, T, depends on the simulation temper- 
ature and the system size and grows dramatically close 
to the critical temperature, Tc, a phenomenon referred 
to as critical slowing down. At Tc the auto-correlation 
time displays a power-law dependence of the system size, 
T ^ L^'^"^ , defining a Monte Carlo dynamical exponent 
zmc- For the well known Metropolis algorithm one esti- 
mate£^'^° zmc ^ 2.1 — 2.2. If efficient algorithms, with 
a very small zmc, cannot be found, this scaling renders 
the Monte Carlo method essentially useless for large lat- 
tice sizes since all data will be correlated. It is therefore 
crucial to develop algorithms with a very small or zero 
Zmc- In order to test the proposed algorithms we shall 
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therefore only consider simulations at the critical point 
since this is where the critical slowing down is the most 
pronounced and where zmc is defined. 

The geometrical worm algorithm.^'^ has proven to be 
very efficient for the study of the bosonic Hubbard model 
when formulated in terms of currents on the links of the 
space-time lattice. The estimated zmc is very close to 
zero*. It is therefore natural to generalize this algorithm 
to classical models defined in terms of bond-variables. 
This is done in section IIIII where a very efRcient algo- 
rithm formulated directly on the dual lattice is devel- 
oped. A related dual algorithm was initially described 
by Prokof 'ev and Svistuno"\ii. Our algorithm can be gen- 
eralized to Potts, clock and other discrete lattice models 
of which the Ising model is the simplest example. The 
dual algorithm also allows for a very simple and elegant 
way of calculating the domain wall free energy directly. 
In section Hill we outline how this is done. Regrettably, 
only in two dimensions is it possible to define such a dual 
algorithm. In section ITTll we also present a directed ver- 
sion of this algorithm where the probability for the worms 
to erase themselves is minimized. However, before pre- 
senting the dual algorithm it is instructive to consider 
geometrical worm algorithms defined on the direct lat- 
tice. Such an algorithm is described in section m in 
both directed and undirected versions. This algorithm 
is applicable in any dimension but is of less interest due 
to its poor efficiency. However, from an algorithmic per- 
spective this algorithm is interesting in its own right and 
it leads naturally to the definition of the dual algorithm. 
Finally, in section Hvl we conclude with a number of ob- 
servations concerning the properties of the algorithms. 



II. GEOMETRICAL WORM ALGORITHMS ON 
THE DIRECT LATTICE 

The first algorithm we present we shall refer to as 
the linear worm algorithm. This algorithm is closely re- 
lated to the geometrical worm algorithms2i2,, however, 
the worms do not form closed loops, instead they form 
linear strings (a worm) of flipped spins. A major advan- 
tage of this algorithm is that we can select the length of 
the linear worm or even the entire distribution of worm 
lengths. This could be advantageous for the study of 
frustrated or disordered models where other cluster algo- 
rithms fail due to the fact that too "big" or too "small" 
clusters are being generated, leading to a significant loss 
of efficiency. 



A. Algorithm A (Linear Worms) 

We begin with a number of useful definitions: In order 
to define a working algorithm we consider two configu- 
rations of the spins, /.t and v related by the introduction 
of a worm. Let /i denote the configurations of the spins 
without the worm, w, and v the configuration of the spins 



with the worm. Furthermore, let Si . . . Sw be the sites 
(or spins) visited by the worm (of length W) and let Ei 
denote the energy require to flip the spin at site Si from 
its position in the configuration /i, with —Ei the energy 
required to flip spin i from its position in v to that in 
ji. Note that Ei is defined relative to the spin configu- 
ration /i. In the following we shall make extensive use 
of the activation or local weight for overturning a spin, 
Af^. More precisely, Af^ denotes the weight for flip- 
ping Si from its position in /i to that in v and A~^^ the 
weight for going in the opposite direction. As we shall 
see, Af^ is not uniquely determined. Here, we shall use 
Af^ = min(l, exp(— AiSi/fcsT)) although other choices 
would be equally suitable. When the worm is moving 
through the lattice it will move from the current site Si 
to a set of neighboring sites cr and it becomes necessary 
to define the normalization Ns^ = Ylia +ct • This nor- 
malization is used for choosing the next neighbor to visit 
among the set cr of neighbors. 

If one considers the proofs for the geometrical worm al- 
gorithmsL^, it is not difficult to see that a generalization 
to classical statistical models defined on the direct lattice 
will depend crucially on the fact that the normalization 
does not depend on the position (up or down) of the 
spin, Si, at the site itself. For the Ising model, defined 
in Eq. J^l, only nearest neighbor interactions are taken 
into account and we can satisfy this requirement simply 
by not allowing the worm to move to any of the 4 near- 
est neighbor sites to the current site Si. In principle, one 
can consider moving the worm to any other sites in the 
lattice, an aspect of this algorithm of considerable intest. 
For the case of longer range interactions we would have to 
restrict the worm to move only to sites that the current 
spin does not interact with. Note that, since we allow 
sites on the same sub-lattice to be visited then the A^^ 
will depend on the order that we visit the sites since sites 
neighboring Si could have been visited previously by the 
worm. As an example of a simple choice for the neigh- 
bors we show in Fig ^ an example where the worm can 
move to 4 neighbors around the current site, excluding 
the 4 nearest neighbors to the site. When defining the 
set of neighbors cr for the worm to move to from Si one 
also has to satisfy the trivial property that if s^+i is a 
neighbor of Si then Si should also be a neighbor of s^+i. 

We now propose one possible implementation of the 
geometrical worm algorithm for the Ising model. A set 
of suitable neighbors to Si is defined as outlined above 
and at the beginning of the construction of the worm we 
draw its length, W from a normalized distribution. We 
refer to this algorithm as linear worms (A). Below we 
outline the algorithm in pseudo-code. 

1. Choose a random starting site, Si=i and with nor- 
malized probability a length W . 

2. Calculate the probability for flipping 
the spin on that site Af^ with Af' = 
min(l,exp(-A£;,,/fcBr)). 
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3. Flip the selected spin, Si. 

4. For each of the k neighbors in the set of neigh- 
bors (T, calculate the weight for flipping, ^f^^ = 
mm{l,exp{-'AEs.^a-/kBT)). Calculate the nor- 
malization Ns- = X](T +0- the probabilities 

5. According to the probabilities pj. select a new site 
Si+i among the neighbors to go to and increase i 
by one, i i + 1. 

6. If i < go toEl 

7. Calculate ^J^"' and Ns„ after the worm 
is constructed and with probability Pdw) — 
1 - mm{l,Af^^NsJ{Aj^«'N,^)) erase the con- 
structed worm. Here Af_^ and -/Vsi are calculated 
before the worm is constructed. 

8. GotoH 

Note that if we decide to construct a worm of length 
= 1 at a given site si then the above algorithm corre- 
sponds to attempting a Metropolis spin-flip at that site. 
See appendix El for a proof of the above algorithm. 

B. Algorithm B (Directed Linear Worms) 

It is an obvious advantage to have control over the 
distribution of the length of the worms. However, if we 
choose the length of the worm at the start of the con- 
struction of the worm as in algorithm A then we allow 
for the worm to backtrack, thereby erasing itself. We 
can try to eliminate or rather minimize the probability 
for the worm to do back-tracking by constructing a di- 
rected algorithm. We closely follow the method outlined 
in reference m in order to construct a directed algorithm. 
See also reference [h| for previous work on directed al- 
gorithms. Let arn and a„ be among the neighbors, a, 
of the site Si. The above proof of algorithm A does not 
depend directly on the definition of the probabilities pj™ 
and Pg" , but only on their ratio, since they have to satisfy 
the following relation: 

This leaves us considerable freedom since we can define 
conditional probabilities, Ps;(m|n), corresponding to the 
probability to continue in the direction at the site Si 
if we are coming from (T„. At a given site we then only 
need to satisfy 

Vsi{m\n) ^ As^+a^ 
Ps.{n\m) As.+cr^' 

If we have k neighbors this defines a fc x fc matrix P where 
the diagonal elements corresponds to the back-tracking 
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FIG. 1: An example of the configuration of the neighbors for 
algorithm A,B. The current site, Si, is indicated by •, the 
neighbors by o. The activation (weight) at the current site 
will depend on the spins at the nearest neighbor sites (shaded 
circles) . 



probabilities, the probabilities for the worm to erase it- 
self. Previously, we had effectively been using ps;(m|n)'s 
which were independent of the incoming direction, (T„, 
since we simply had ps^ (rn|n) = Ag^+o-^ 1^ Si , correspond- 
ing to a matrix P with identical columns. However, the 
above condition of balance, Eq. ^ leaves sufficient room 
for choosing very small back-tracking probabilities and in 
most situations the corresponding diagonal elements of P 
can be chosen to be zero. If we impose the constraint that 
the sum of the diagonal elements of P should be mini- 
mal the problem of finding an optimal matrix P can be 
formulated as a standard linear programming problem to 
which conventional techniques can be applied**. 

At a given site Si we choose 4 neighbors, as indicated 
in Fig. n by the open circles o. The activation at a given 
neighbor will depend on its 4 nearest neighbors shown as 
the shaded circles in Fig. ^ Technically, each p3;(TO|n) 
now depends on the position of all the 16 surrounding 
spins, and hence there are in principle 2^® possible ma- 
trices P at each site. It is therefore not feasible to choose 
too large a set of neighbors when directing the algorithm. 
In the absence of disorder, the 2^^ matrices can be tab- 
ulated and minimized at the beginning of the simulation 
and for the isotropic ferromagnetic model at hand it is 
easy to see that at most 625 different matrices occur. 

We can now use the above matrices P for an algorithm 
that performs directed linear worms of length varying be- 
tween 1 and W . We again assume that a set cr of A: 
neighbors has been chosen. 

1. Choose a random starting site, Si=\ and with nor- 
malized probability a length W . 

2. Calculate the probability for flipping 
the spin on that site A^^ with A^^ — 
mm{l, eM~^EsJkBT)). 



4 



3. Flip the selected spin, Si. 

4. If Si — si then for each of the k neighbors in 
the set (T, calculate the probability for flipping, 



A 



i+t7 



= min(l, exp(— Ai?s.-|_cr/fcsT)). Calculate 



the normalization Ns- = ^fi+a the prob- 
abilities Pg. = ^f'i|_£,/A^si • Else: Depending on 
the configuration of the nearest neighbors of the k 
neighbors select the correct (minimized) matrix P 
and if the worm arrived from direction (t„ set p^. 
equal to the n'th column of P. 

5. According to the probabilities pj. select a new site 
Si-f 1 to go to and increase i by one, i — > « -I- 1. 

6. If i<W go toEl 

7. Calculate Ar^"' and 7V,,„ after the worm 
is constructed and with probability Pe{w) — 
1 - mm{l,Af^^NsJiAj^^Ns^-)) erase the con- 
structed worm. Here Af_^ and Ns^ are calculated 
before the worm is constructed. 

8. GotoH] 

The proof of the above algorithm can be found in ap- 
pendix 



C. Performance: Algorithms A,B 
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FIG. 2; Autocorrelation functions for the energy Csit) as 
a function of Monte Carlo time, calculated at Tc- Shown 
are results for L — 10 for the Metropolis, Swendsen-Wang 
and for the linear worm (A) and directed linear worm (B) 
algorithms (with uniform worm length). Note that in order 
to compare the 4 algorithms the time axis has been scaled 
so that in all cases 1 time step corresponds to an attempted 
update of all the variables. The autocorrelation functions 
shown were calculated averaging over 20 million worms (MCS 
for the metropolis algorithm). 



In order to test the performance of algorithms A and 
B we consider the autocorrelation function Co{t) of an 
observable O. This function is defined in the standard 
way: 



Co{t) 



{o{t)om - {Of 

(O^) {Of 



(4) 



For a reliable estimate of Co {t) we typically generate 20 
million worms. In Fig. [21 we show results for the auto- 
correlation function for the energy C'sit) for the linear 
worm algorithms A, B for a system of size L = 10. In both 
cases the worm length W was chosen from a uniform dis- 
tribution. This is compared to Csit) for the usual single 
flip Metropolis algorithm and the Swendsen-Wang algo- 
rithm. For the latter two algorithms the Monte Carlo 
time is usually measured in terms of Monte Carlo Sweeps 
(MCS) where an attempt to update all the spins in 
the lattice has been made. However, for algorithms A 
and B a worm will on average only attempt to update 
(W) of the spins. Hence, if time is measured in terms 
of generated worms for the worm algorithms it should 
be rescaled by a factor of {W)/L'^ for a fair comparison 
to be made with algorithms where Monte Carlo time is 
measured in MCS. This has been done in Fig. |21 From 
the results in Fig. |2it is clear that the efficiency of algo- 
rithm A and B is fairly poor and worse than the much 
simpler Metropolis algorithm. We have checked that this 
remains true for significantly larger lattices sizes. 



The main cause of this poor behavior is the necessity 
to include a rejection probability Pg. If a worm on aver- 
age attempts to update (W) spins, we need on average to 
generate L'^/{W) worms to complete a MCS. As a rough 
estimate, let us consider that the worm is accepted with 
the average probability p and that only a fraction x of the 
(W) attempted spin flips result in an updated spin (a spin 
can be visited several times) . It then follows that on av- 
erage pxL'^ are actually flipped in a MCS with the linear 
worm algorithm. If the average probability for flipping 
a spin with the Metropolis algorithm is q and if all the 
spins selected in a MCS are different, the corresponding 
number for the Metropolis algorithm is qL^, presumably 
of the same order as pxL^ and likely larger. The effect 
of the directed linear worm algorithm is to maximize x 
as much as possible and it therefore seems unlikely that 
the linear worm algorithm in its present form can per- 
form better than the Metropolis algorithm unless p also 
is maximized. 

It would be very interesting if algorithms A and B 
could be modified so that p = 1 {Pe = 0). We have so far 
been unable to do so. Such a modified algorithm would 
be significantly more efficient than the Metropolis algo- 
rithm (but presumably less efficient than the Swendsen- 
Wang algorithm). It would be much more versatile and 
could be of significant interest for frustrated or disordered 
systems since it would not require the construction of 
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clusters but only of the much simpler worms, the length 
of which can be chosen. 

Even though it is quite interesting to be able to choose 
the distribution of the worm lengths at the outset, the 
question arises which distribution of worm lengths will 
give the most optimal algorithm. In the present work we 
have chosen a distribution of worm lengths that is uni- 
form between 1 and 2L, but on general grounds we expect 
a power-law form for this distribution to be more opti- 
mal at the critical point. For the study of Bose-Hubbard 
models using geometrical worm algorithms it was noted^ 
that the distribution of the size of the worms follow a 
power-law with an exponent of approximately 1.37 at the 
critical point. Hence, in the present case it would appear 
likely that an optimal power-law distribution of the worm 
lengths at Tc can be found, defining a new "dynamical" 
exponent. The idea of choosing the distribution of the 
worm size to optimize the algorithm resembles previous 
work by Barkema and Newmani^ti^. 

In closing this section, we note that it is quite straight- 
forward to define an algorithm on the direct lattice where 
the worms form closed loops (rings). In this case the 
length of the worm is determined by the size of the loop. 
We have tested such algorithms but their performance is 
even worse than algorithms A and B since the constraint 
that the initial spin has to be revisited makes the length 
of the worms in some cases diverge, or for other variants 
of such an algorithm, go to zero. 

III. DUAL ALGORITHM 

It is clear that one problem with both of the above 
algorithms is the fact that the spin on the initial site is 
treated differently than the remaining spins. This is be- 
cause the geometrical worm algorithms are more suitable 
for an implementation directly on the dual model. Hence, 
we will now try to describe an algorithm that moves a 
worm along the dual lattice by updating bond-variables 
on the direct lattice. 

We begin with some definitions analogous to the treat- 
ment of KadanofEi^. On the direct lattice we define, at 
each site (j, k) an integer variable Sj,k = ±1. Here j de- 
scribes the index in the x-direction and k the index in 
the y-direction. Then we can define the following bond- 
variables: 

k) = Sj+i.kSj.k (5) 
b^{j,k) = Sj.k+iSj^k- (6) 

The Ising model has variables where as we see that we 
have 2L^ bond- variables. However, it is easy to see that 
the bond-variables satisfy a divergence free constraint at 
each site of the dual lattice: 

b'^ij, k) + + 1, fc) - h'^ij, k + 1)- b^U, A:)(mod4) = 

(7) 

giving us constraints. However, if we define the model 
on a torus the constraint on each dual lattice site is equal 



to the sum over the constraints on all the other dual 
sites. Hence, in this case we obtain only — 1 indepen- 
dent constraints. In addition we also have to satisfy the 
boundary conditions, SL+i,k = si^k, = Sj,!- This 

implies that for an L x L lattice with even L: 

L 

k)-L (mod4) = V/c (8) 

L 

^6''(j,fc) -L (mod4) = OVj. (9) 

k=l 

These constraints are not independent of the previously 
defined constraints Eq. j?)). In fact, it is easy to see 
that if the boundary constraints Eq. are applied at 
just one row and one column then the divergence free 
constraints Eq. ^ will enforce the boundary constraints 
at the remaining rows and columns. Hence, these con- 
straints only give us 2 more independent constraints, in 
total L'^ + 1 constraints. The model, written in terms 
of the bond- variables, therefore has — 1 free variables 
and correspondingly half the number of degrees of free- 
dom compared to the formulation in terms of the spin- 
variables. This is a natural consequence of the fact that 
the bond-variable model does not distinguish between a 
state and the same state with all the spins reversed. 

The partition function for the Ising model on a L x L 
torus can now be written in the following manner: 

Z = Trfc.Trfc„ 11^=1 IlLi 
X exp {K-kb^ij, k) + Kl^by{j, fc)} . (10) 

Here Tr' denotes the trace over bond-variables satisfy- 
ing the above constraints and f. = J^^/ksT, = 

A. Algorithm C (Dual Worm Algorithm) 

We now turn to a discussion of the dual algorithm. 
From the above description in terms of bond- variables it 
is now quite easy to define a geometrical worm algorithm 
on the dual lattice closely following previous work on such 
algorithmsSi^. We denote the i'th site on the dual lattice 
that the worm visits as di. A "worm" is constructed 
by going through a sequence of neighboring sites on the 
dual lattice by each time choosing a direction a to follow 
according to an appropriately determined probability. 
When the worm moves from di to c?i+i the corresponding 
bond- variable bi that the worm crosses is flipped. The 
bond-variables can take on only 2 values ±1. Hence, 
the associated energy cost for flipping hi is given by 
AE — 2 J"j,6"(j, fc), a = x,y. We can now define weights 
for each direction cr, A^" — min(l, exp(— Ai^'^/fesT)) 
used for determining the correct probability for choosing 
a new site. This is not the only choice for the weights. 
Other equivalent choices should work equally well. The 
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FIG. 3: A worm moving through the dual lattice. The direct 
lattice is indicated by solid lines and the spins on the direct 
lattice by solid circles. The dual lattice is indicated by dashed 
lines and the sites on the dual lattice by open circles. As the 
worm moves through the lattice the bond- variables are flipped 
as indicated by the thick bonds in the figure. 
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FIG. 4: An illegal worm move. The worm winds around the 
lattice in the vertical direction and correspondingly all rows 
have an odd number of flipped bond-variables 6. The direct 
lattice is indicated by solid lines and the spins on the direct 
lattice by solid circles. The dual lattice is indicated by dashed 
lines and the sites on the dual lattice by open circles. As the 
worm moves through the lattice the bond- variables are flipped 
as indicated by the thick bonds in the flgure. 



bond-variables are updated during the construction of 
the worm and the worm is finished when the starting site 
on the dual site is reached again and the worm forms 
a closed loop. The dual worm algorithm can then be 
summarized using the following pseudo-code: 

1. Choose a random initial site di on the dual lattice. 

2. For each of the directions a = ±x,zty calculate 
the weights A-^" associated with flipping the bond- 
variable perpendicular to that direction, A^" =^ 
min(l, exp(-A£;'"/A:BT)). 

3. Calculate the normalization Nd, = J2a ^^'^ 
associated probabilities p^. ~ A^" jN^^- 

4. According to the probabilities, , choose a direc- 
tion a . 

5. Update the bond- variable hi for the direction cho- 
sen and move the worm to the new dual lattice site 

6. If di ^ rfi gotoEl 

7. Calculate the normalizations Ndx a-nd Ndx of the 
initial site, si, with and without the worm present. 
If the worm is "legal" , i.e. with even winding num- 
ber in both the x- and the y-direction (O^, 
both - see definition below), then erase the worm 
with probability Pg = 1 - min(l, iVd^/TVdJ. If the 
worm is "illegal", that is if cither or cal- 
culated when the worm has closed equal 1, then 
always erase it. GotoQ] 

Following the above discussion of the boundary con- 
straints it is easy determine if the winding number in the 
X- or y-direction is odd by simply choosing one row /cq 
and one column jo and calculating the number of frus- 
trated, 6 = — 1, bond- variables: 

= ^ l^fe^^lM (,,od2) 

i 

0^-E^^^(-d2), (11) 

since a worm with an odd winding number in the y- 
direction will result in being 1 independent of fco, 
or equivalently for the x-direction and . See Figure 01 
and then takes on the values 1 or depending on 
whether the corresponding winding numbers are odd or 
even. 

Several points are noteworthy about this algorithm. 
To a great extent it is simply the dual version of the 
Wolff algorithm^. Each worm encloses a cluster of 
spins on the direct lattice (not necessarily a simply 
connected cluster). Flipping the bond- variables in the 
worm effectively flips the spin in the cluster. This corre- 
spondence is particular to two dimensions which is the 
only dimension where the present dual worm algorithm 
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can be defined. The intuitive argument for this is that 
the bond-variables updated by the worm have to enclose 
a finite volume of spins on the direct lattice, this is 
only possible in two dimensions. In dimensions higher 
than two the one dimensional worms are not capable of 
enclosing a finite volume. Mathematically, it can be seen 
that the divergence less constraint, Eq. (71), cannot be 
satisfied under worm moves. Under most conditions the 
rejection probability, Pe, is very small, however, for small 
lattice sizes the probability for generating a worm with 
an odd winding number in either the x- or y-direction 
can be non-negligible. It is important to note that since 
the bond-variables are updated during the construction 
of the worm the generated configurations are not valid 
and the above constraints are not satisfied. However, 
once the construction of the worm is finished and the 
path of the worm closed, the divergenceless constraint is 
satisfied. Rejecting worms with odd winding numbers 
then assures that the resulting configurations are valid. 
Also, when the worm moves through the lattice it may 
pass many times through the same link and cross itself 
before it reaches the initial site where the construction 
terminates. Finally, at each step i in the construction 
of the worm it is likely that the worm at the site di will 
partially "erase" itself by choosing to go back to the site 
di-i visited immediately before, thereby "bouncing" off 
the site di. 

Proof of Algorithm C (Dual Worm Algorithm) 

Now we turn to the proof of detailed balance for the al- 
gorithm. As before, let fj, denote the configuration of the 
bond-variables without the worm, w, and v the configu- 
ration with the worm. Let us consider the case where the 
worm, w, visits the sites {di . . . dw} on the dual lattice. 
Here di is the initial site. The worm then goes through 
the corresponding bond-variables {61 . . . bw}- Note that, 
dw is the last site visited before the worm reaches c?i. 
Furthermore, let Ei denote the energy required to flip the 
bond-variable bi from its position in the configuration /i, 
with —Ei the energy required to flip bi from its position 
in to that in /i. The total probability for constructing 
the worm w is then given by: 



P{w; di dw) 



w 



P(di)(l-PeM)P(legal|w;)n 



(12) 



The index a denotes the direction needed to go from di 
to c?i+i, (ice, ±y), crossing the bond-variable bi. -P(di) is 
the probability for choosing site di as the starting point 
and Pe{'w) is the probability for erasing a "legal" worm 
after construction. Finally, P(legal|w) is the conditional 
probability that the constructed worm w is "legal" , with 
even winding numbers in both x- and y-directions. If 
the worm w is legal and has been accepted we have to 
consider the probability for reversing the move. That 
is, we consider the probability for constructing an anti- 



A d3 



a 

A 



d2 



h 

bw 

■_--^a-- 

bl d/ \n 



I 

o- 

A 



,dw-2 



bw-l 



l^bw-2 



6 



\ ^dw-1 



FIG. 5: The notation used for the proof of the dual worm algo- 
rithm. The configuration shown corresponds to a worm (solid 
line) partly erased by its corresponding anti-worm (dashed 
line) . 



worm w annihilating the worm w. Note that, in this case 
the sites are visited in the opposite order, d\ — d\^di = 
dw,...,dw = d2, in general di = dw-i+2 (* 7^ 1)- In 
the same manner the bond-variables are also visited in 
the reverse order bi = bw,b2 = bw-i, ■ • • ,bw = bi, in 
general bi = dw-i+i- Consequently, when describing 
w, we directly use d, b instead of d, b. The notation is 
illustrated in Fig. ]E\ We therefore have: 

Piw;di^dw) = P{di){l - Peiw))P{legal\w) 



N, 



n 



1 i=W 



(13) 



Please note that the index i runs over W,...,2. Let us 
now consider the case where both of the worms w and w 
have reached the site di different from the starting site 
di . See Fig. |S1 Since we are updating the link variables 
during the construction of the worm we immediately see 
that = Na, i I. Furthermore, A^* and A'^' 
only depend on the bond- variable bi, and we therefore 
see that: 



A' 



A 



1 . . . I^. 



(14) 



Hence, since P{di) is uniform through out the dual lat- 
tice, and since w and w must wind the lattice in precisely 
the same manner resulting in P(legal|w) — P(legal|w), 
we find: 



PH ^ 1 - P.jw) Ng, 
Piw) l-P,{w)Nd, 



eM-^ETot /keT), (15) 
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where Ai?Tot is the total energy difference between a con- 
figuration with and without the worm w present. With 
our definition of Pe we see that (1 — Pe(w))/(1 — Pe(w)) = 
I ■ Hence, with this choice of Pe we satisfy detailed 
balance since: 

^ = expi^AETot /ksT). (16) 

As was the case for the previous algorithms there are 2 
ways of introducing a worm that will take us from fi — s- 
I/, and equivalently 2 ways of introducing ui. (Strictly 
speaking, each of the two ways actually represents an 
infinite sum as discussed under the proof of algorithm A) . 
In both cases Eq. IjKil) holds and it follows that P(/i 
v) = P{v /i) cxp(— Ai?Tot/fcsr) showing that detailed 
balance is satisfied. Ergodicity is simply proven as the 
worm can perform local loops and wind around the lattice 
in any direction. ■ 

Illegal worms: One should be cautious when treat- 
ing "illegal" moves such as the worms with odd winding 
numbers encountered above. Suppose that our Monte 
Carlo algorithm has generated a number of configura- 
tions Ci, C2, C2, C3, . . . , Ci_i, Ci where some configura- 
tions are repeated since the worm moved has been re- 
jected. From the configuration Ci we now construct a 
worm that turns out to be "illegal" with odd winding 
numbers. From the above proof we then see that we have 
to repeat Ci in the sequence of configurations since not 
repeating the configuration would lead to a total accep- 
tance probability different from (1 — Pe(w))P(legal|w). 
For instance, imagine we kept generating worms till we 
found a "legal" one. This would make the total accep- 
tance ratio a sum over the number of times we try which 
would be incorrect. 



B. Algorithm D (Dual Directed Worm Algorithm) 

Following the discussion of algorithm B it is straight 
forward to develop a directed version of algorithm C. As 
was the case for algorithm B we define conditional prob- 
abilities P(i;(TO|n), corresponding to the probability for 
continuing in the direction cjm if the worm has come from 
direction (t„. Since we in algorithm C, at each site on the 
dual lattice d^, can go in four directions this leads us to 
define a 4 x 4 matrix P at each site on the dual lattice. 
However, compared to algorithm B, the number of differ- 
ent matrices at a given dual site is now much smaller since 
the weights A^^ only depend on the associated bond- 
variable. In fact there are now only 16 possible matrices 
at a given site and these can easily be optimized and 
tabulated at the outset of the calculation. Even in the 
presence of disorder, a large number of them can be tabu- 
lated. This directed dual worm algorithm is now identical 
to algorithm C except for the fact that if di is different 
from o?i, the pJJ. are selected from these optimized matri- 
ces (the same way at it was done for algorithm B). Due 
to the simplicity of this modification and the similarity 




t <W>/(2 L^) 



FIG. 6: Autocorrelation functions for the energy Csit) as a 
function of Monte Carlo time, calculated at Tc. Shown are 
results for L = 10, 20 for the Metropolis, dual worm and dual 
directed worm algorithms. Note that, in order to compare 
the 3 algorithms the time axis has been scaled so that in all 
cases 1 time step corresponds to an attempted update of all 
the variables. The inset shows the scaling of the autocorre- 
lation time r' as a function of the system size for the dual 
worm and dual directed worm algorithm. The autocorrela- 
tion time is here rescaled: r' = r(W)/(2I/^). In both cases 
we find roughly a power-law form r' ~ L'''^® shown as the 
solid lines in the inset. The autocorrelation functions shown 
were calculated averaging over 10-20 Million worms (MCS for 
the metropolis algorithm). 



with algorithm B, we do not present pseudo-code nor a 
proof for this directed algorithm. It is easy to see that 
rejection probability remains identical to the one used in 
algorithm C, Pe{w) = 1 - inm{l,NdjNd,)- 



C. Performance Algorithm C,D 

In order to compare the two algorithms C and D de- 
fined on the dual lattice we have again calculated the 
autocorrelation function CE(t). Our results are shown 
in Fig. and tabulated in Table HJ As on the direct 
lattice we use the energy to calculate the autocorrela- 
tion functions. For algorithms C and D we have con- 
structed 10 million worms for the smaller lattices and 
20 million for the largest lattice (L=160). We compare 
to the Metropolis algorithm where each time step cor- 
responds to an attempted update of the spins. Since 
the worm algorithm in this case on average only attempts 
to update (W) out of the 2L^ bond-variables we have 
to rescale the time-axis by a factor of (VF)/(2L^) when 
we compare algorithms C and D to the Metropolis algo- 
rithm. We can therefore define an autocorrelation time r 
in units of worms constructed and a rescaled autocorre- 
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TABLE I: Average worm size defined as the average num- 
ber of bond-variables attempted to change during the con- 
struction of a any worm, including worms eventually rejected. 
Also listed is the estimated autocorrelation time, t, in units 
of worms constructed and the rescaled autocorrelation time 
r' = r(W>/(2L^) in units of 2L^ attempted updates. The 
autocorrelation time is here here defined by Ce(r) = 0.05. 
Results are listed for the dual and dual directed worm algo- 
rithms. 
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29.3 
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25.7 
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30 
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31.7 
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66.6 


14.8 


40 
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37.0 
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81.9 


16.9 


80 
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303.8 


50.3 


2204 


135.9 


23.4 


160 


7132 


478.8 


66.7 


7408 


214.3 


31.0 



lation time r' = r(W)/(2L^) in units of attempted 
updates. In order to simplify the analysis of the data 
we define the autocorrelation time as Ce{t) = 0.05. As 
can be seen in Fig. |H| the dual worm algorithms repre- 
sents a dramatic improvement over the Metropolis algo- 
rithm. We have not plotted results for the Swendsen- 
Wang algorithm on the direct lattice since for L — IQ 
they are indistinguishable from the results obtained for 
the dual directed worm algorithm. The scaling of r' with 
L is a relatively small power-law r' ^ LP '^^ . This expo- 
nent is small but larger than most values quoted for the 
Swendsen-Wang and Wolff algorithms. This is perhaps 
not too surprising since the worms are allowed to cross 
themselves thereby erasing previous updates. The differ- 
ence in autocorrelation times between algorithm C and 
its directed version D appears to be a constant factor. 
Therefore, the direction of algorithm C does not change 
the exponent zmc- 

For a 10 X 10 lattice at Tc roughly 85% of all con- 
structed worms are accepted and about 10% of the con- 
structed worms are "illegal" when simulations are per- 
formed using the dual geometric worm algorithm. For 
the directed version of the algorithm the corresponding 
numbers are 75% and 20%. The number of illegal worms 
drops rapidly with L and for L = 30 90% of the worms 
are accepted with 5% illegal worms for the undirected al- 
gorithm compared to 83% and 10% for the directed dual 
algorithm. 

The presented dual algorithms are quite general since 
they depend only in a relatively limited way on the under- 
lying Hamiltonian. Even though they may not be com- 
petitive with the Swendsen-Wang and Wolff algorithms 
for the Ising model they may be quite interesting to apply 
to other two-dimensional models that can be formulated 
in terms of bond- variables and where more effective clus- 
ter algorithms are difficult to define. 



D. Domain Wall Free Energy 

A quantity of particular interest in many studies of 
ordering transitions is the domain wall free energy, the 
difference in free energy between configurations of the 
system with periodic (P) and anti-periodic (AP) bound- 
ary conditions in one direction, AF = Fap — i^p. The 
domain wall free energy can be shown to obey simple scal- 
ing relations and is a very efficient tool for distinguishing 
between different ordered phases. Unfortunately, it is 
usually very difficult to calculate free energies directly 
in Monte Carlo simulations. A clever trick to do so 
is the boundary flip MC method, proposed by Hasen- 
buscbiSiik, where the coupling constants at the boundary 
are considered as dynamical variables, and can be flipped 
during the course of the MC simulation. If Pp{T) and 
PaP{T) describe the probability to obtain MC configu- 
rations with periodic and anti-periodic boundary condi- 
tions and and Zp, Zap the associated partition functions, 
then the domain wall free-energy is given by: 

g-/3(FAP-Fp) ^ ^AP ^ Pap{T) ^^^^ 

Zp Pp (T) 

Effectively, the boundary flip MC method attempts to 
sample the ratio Z^p/Zp^ and this is usually cumber- 
some. However, for the dual algorithms it is trivial to 
obtain this ratio. It follows from Eq. 1^, that if we al- 
low the constructed worms to have all different kinds of 
winding around the torus (making all worms legal) then 
we are sampling a partition function which is a sum of 
four terms: 

^ — ^Px,Py + ^APx.Py + ^Px.APy + ^AP3j,APy (18) 

Here, ZBCx,BCy is the partition function with BC2:(BCy) 
in the x(y)-direction. This is so, because we can turn an 
illegal worm into a legal worm by changing the bound- 
ary condition in the direction(s) where the winding of 
the worm violates Eq. @. During the MC simulation 
it is easy to see to which term a given configuration 
contributes simply by keeping track of the total wind- 
ing number in both the x- and y-direction. See eq. Hll|l . 
When this is even, the boundary condition in the associ- 
ated direction is periodic, and if it is odd the bound- 
ary condition is anti-periodic. In order to calculate 
■Z^APx.Py/'Z^Px.Py ^6 therefore simply have to count how 
many times the total winding number is odd in the x- 
direction and even in the y-direction and divide with the 
number of times it is even in both direction. Due to its 
simplicity, we expect this approach to be an extremely 
efficient way of calculating AF. In tabled we show pre- 
liminary results for AF/kT calculated with 10^ worms 
for the two-dimensional Ising model at Tc- For reference 
we compare to exact results obtained using pfaffiansii. 

We are grateful to Philippe de Forcrand for suggesting 
the above way of calculating AF. 
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TABLE II: The domain wall free-energy AF/kT for several 
system sizes calculated at the critical temperature Tc obtained 
using our MC method and exact results from Ref. ilTl 



_L Exact MC 

~4 0.9658246670 0.9654(5) 

16 0.9853282229 0.985(10) 

24 0.9859725679 0.985(10) 

32 0.9861972559 0.986(10) 

48 0.9863574947 0.986(10) 

64 0.9864135295 0.984(20) 



to a glassy phase). One might have expected the hn- 
ear worm algorithm to perform significantly better than 
the Metropolis algorithm when disorder is present. We 
assume that the fact that it is only comparable to the 
Metropolis algorithm is because large parts of the spins 
are not effectively flipped since worms are mostly rejected 
in those parts of the lattice. If this is true, a modified 
version of the linear worm algorithm with Pg = or with 
a much more ingenious choice of the neighbors would be 
of considerable interest since it would hold the poten- 
tial to sample such disordered models significantly more 
effectively than the Metropolis algorithm. 



IV. DISCUSSION 

We have presented two generalizations of the geometri- 
cal worm algorithm applicable to classical statistical me- 
chanics model. The first algorithm exploits linear worms 
on the direct lattice, the other closed loops of worms on 
the dual lattice. In both cases have we developed directed 
versions of the algorithms. While the first algorithm is 
applicable in any spatial dimension, it is quite clear that 
for topological reasons it is only possible to define a worm 
algorithm of the proposed kind on the dual lattice in two 
spatial dimensions. 

Due to its poor performance the linear worm algorithm 
(A) is mainly of interest from the perspective of fostering 
the development of new and more advanced algorithms. 
The linear worm algorithm presents a number of very 
attractive features such as the possibility to choose the 
distribution of the worm lengths and it seems possible to 
significantly improve the performance of the algorithm 
(and most other geometrical worm algorithms) by elim- 
inating the rejection probability. So far, we have been 
unable to do so, but this would clearly be very interesting 
since the algorithms then would be quite promising for 
the study of frustrated or even disordered models where 
other cluster algorithms fail. 

The dual worm algorithm is very efficient and even 
though the obtained zmc is larger than for the Swendsen- 
Wang algorithm it is quite competitive since it requires 
relatively little overhead and is very simple to implement. 
One of the most interesting features of this algorithm is 
the topological aspect of the worms, resulting in forbid- 
den worm moves. 

The geometrical worm algorithm has in some cases 
been very successfully applied to the study of models 
with disorder—. We therefore applied both the linear and 
dual worm algorithm to the bond disordered Ising model. 
In both cases do the algorithms perform worse than the 
Metropolis algorithm. For the dual worm algorithm this 
failure is due to the fact that the average size of worms 
diverge when glassy phases are approached. The linear 
worm algorithm would seem more promising for study 
of models with disorder since the length of the worms 
can be controlled from the outset. In fact the perfor- 
mance of the algorithm is in this case quite comparable to 
the Metropolis algorithm (which has a diverging t close 
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APPENDIX A: PROOF OF ALGORITHM A 

The probability for creating the worm, w, starting from 
site si is given by: 



-(1-PeM) (Al) 



Here P{W) denotes the probability for choosing a worm 
length of W . The probability for introducing an anti- 
worm, starting at the site si, erasing the worm by 
going in the opposite direction along w becomes: 



P{w-sw~^si)= P{W)P{sw) 



a; 



'Ew -1 



-(l-Pe(^)))(A2) 



Since neither Ns- nor Ns- depend on the position of the 
spin on the site i and since the spins are visited in the 
precise opposite order for the anti-worm with respect to 
the worm, we then see that Ng- = Nst Vi 7^ 1, W, and we 
find: 

P{w; SI -> sw) _ 1 - Peiw) ...Af^ 7V,„ 



P{w; sw ^ si) 1 - P,iw) . . . ^-^^1 N,, 

(A3) 

Since Af.'/A^.^^ = exp{— AEi / ksT) we see, using our 
definition of P„, that 



P{w; si Sw) 
P{w;sw -> si) 



exp(-A£;/fcsr). (A4) 
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Quite generally, there is more than one way to introduce 
a worm to go from /j, to v. In fact there are 2 ways of 
introducing a worm taking the system from configuration 
fj, to V, one where the worm traverses the sites from si to 
sw and another where the worm traverses the sites from 
sw to si. Equivalently, there are 2 ways of introducing 
the anti-worm. Hence, v) is, a sum of the two 

probabilities, P(/z ^ v) = P{w; si sw) + Piw; sw 
si). However, employing the result we have just proven 
we see that: 



matrix P, at the site Si for continuing to the site s^+i 
knowing that the worm is coming from the site Si_i. 
Starting from the configuration v we can now calculate 
the probability for introducing an anti-worm, w, starting 
at the site sw, erasing the worm, w. We find: 



P{w; sw^si )^P{W)P{sw) 



A. 



-Psiv-lisW-2\sw) 



Ps3(s2|s4)Ps2(sik3)(l - Peiw)). 



(B2) 



P{^.^v)= [P{w;sw ^ si) + P{w;si ^ Sw)] 
X exp{-AE/kBT) 
= P{i^ ^ ^i)exp{-AE/kBT). (A5) 

Hence, the proposed algorithm satisfies detailed balance 
on a bi-partite lattice and we see that detailed balance 
is satisfied no matter how many possible ways one can 
introduce a worm in order to go from /x — > z/ as long 
as the number of possible anti-worms matches. Strictly 
speaking, since the worm can partly erase itself, there 
is in fact an infinite number of worm-moves that will 
take us from /i to however, in each case there is a 
matching anti-worm so we can replace the above sum 
over two probabilities with an infinite sum, yielding the 
same conclusion. Ergodicity is proven by noting that 
single spin flips are incorporated in the algorithm. ■ 



APPENDIX B: PROOF OF ALGORITHM B 

The proof is quite similar to the proof of algorithm 
A and we use the same notation. The probability for 
creating a worm, w, starting from site si is then given 
by: 



P{w;si 



sw)^pmp{^i)^vsA^zVi)- 



N. 



Si 



Ps^v-2isW-l\sW-3)Psw-lisw\sW-2)-{'^ - PeM^'Bl) 

Here P{si) denotes the probability for choosing site si 
as the starting site for the worm, P{W) the probability 
for choosing a worm length W. (si+i |si„i) denotes 
the conditional probability, extracted from the minimized 



At a given site, Si, the configurations of the spins during 
the construction of the worm and the anti-worm only 
differ by the position of the spin at the site itself. Since 
the set of neighbors a exclude any nearest neighbors to 
Si we see that Pg^ = Psi ■ Hence, by construction. 



Psi{Si+i\Si-i) _ A^V^/ 



(B3) 



Using the definition of and the fact that P(si) 
P{sw), we then see that: 



Piw;si^sw)_ Af^^...Af- 



P{w;sw'^si) A-f'...A7^' 



(B4) 



Since Af.' /Aj.^' = cxp{-AEi/kBT) wc finally see that 



P{w;si Sw) 
P{w; sw si) 



exp(-A£;/fcBr), 



(B5) 



where AE is the total energy cost in going from configu- 
ration fi to configuration v. As was the case for algorithm 
A, there are 2 ways of introducing a worm taking the sys- 
tem from configuration fi to v (modulo sites visited an 
even number of times), one where the worm traverses the 
sites from si to sw and another where the worm traverses 
the sites from sw to si. Equivalently, there are 2 ways 
of introducing the anti-worm. We again see that detailed 
balance is satisfied. Ergodicity is proven by noting that 
single spin flips are incorporated in the algorithm, since 
for W — I the algorithm corresponds to attempting a 
Metropolis spin-fiip at the selected site.B 
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